Visualization of the 2020 Environmental Performance Index (EPI)

Load libraries

library("reshape2")
library("ggplot2")
theme_set(theme_bw(base_size = 12) +
            theme(rect = element_rect(fill = "transparent")))
library("leaflet")
library("rgdal")

## Loading required package: sp

## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.4.0, released 2018/12/14
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
## Path to PROJ shared files: (autodetected)
## Linking to sp version:1.4-5

library("htmlwidgets")
library("lavaan")

## This is lavaan 0.6-7

## lavaan is BETA software! Please report any bugs.

library("lavaanPlot")

Read data

Clear workspace, set working directory

rm(list = ls())
setwd("~/epi")

Load helper functions

These helper functions retrieve, read and process data.

source("UtilityFunctions.R")

Import EPI data

Run the helper function importEPI() and release the contents of the resulting list into the global environment.

epi.list <- importEPI()

## 'data.frame':    8234 obs. of  7 variables:
##  $ code         : int  4 24 8 784 32 51 28 36 40 31 ...
##  $ iso          : chr  "AFG" "AGO" "ALB" "ARE" ...
##  $ country      : Factor w/ 179 levels "Afghanistan",..: 1 4 2 170 6 7 5 8 9 10 ...
##  $ region       : Factor w/ 8 levels "Asia-Pacific",..: 7 8 2 5 6 3 6 4 4 3 ...
##  $ X2019        : num  1.93e+10 8.88e+10 1.53e+10 4.21e+11 4.45e+11 ...
##  $ EPI.new      : Factor w/ 46 levels "AGR","AIR","APE",..: 12 12 12 12 12 12 12 12 12 12 ...
##  $ EPI.new.value: num  25.5 29.7 49 55.6 52.2 52.3 48.5 74.9 79.6 46.5 ...

list2env(epi.list, .GlobalEnv)

## <environment: R_GlobalEnv>

rm(epi.list)

Retrieve world borders

Download the shape file, if it is absent. Note: overwrites existing files!

if(!file.exists("2020/TM_WORLD_BORDERS-0.3.shp")){
  print("File does not exist, downloading it now:")
  download.file("thematicmapping.org/downloads/TM_WORLD_BORDERS-0.3.zip",
                destfile="world_shape_file.zip")
  system("unzip -u world_shape_file.zip;
         rm world_shape_file.zip Readme.txt")
  }

The, read this shape file into a SpatialPolygonsDataFrame.

world.spdf <- readOGR(dsn = "2020", layer = "TM_WORLD_BORDERS-0.3",
                      verbose = TRUE)

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/r/epi/2020", layer: "TM_WORLD_BORDERS-0.3"
## with 246 features
## It has 11 fields
## Integer64 fields read as strings:  POP2005

summary(world.spdf)

## Object of class SpatialPolygonsDataFrame
## Coordinates:
##    min      max
## x -180 180.0000
## y  -90  83.6236
## Is projected: FALSE 
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Data attributes:
##      FIPS               ISO2               ISO3                 UN       
##  Length:246         Length:246         Length:246         Min.   :  4.0  
##  Class :character   Class :character   Class :character   1st Qu.:215.0  
##  Mode  :character   Mode  :character   Mode  :character   Median :429.0  
##                                                           Mean   :431.8  
##                                                           3rd Qu.:650.5  
##                                                           Max.   :894.0  
##      NAME                AREA             POP2005              REGION      
##  Length:246         Min.   :      0.0   Length:246         Min.   :  0.00  
##  Class :character   1st Qu.:     44.5   Class :character   1st Qu.:  2.00  
##  Mode  :character   Median :   5515.5   Mode  :character   Median : 19.00  
##                     Mean   :  52696.1                      Mean   : 65.43  
##                     3rd Qu.:  34708.8                      3rd Qu.:142.00  
##                     Max.   :1638094.0                      Max.   :150.00  
##    SUBREGION           LON               LAT          
##  Min.   :  0.00   Min.   :-178.13   Min.   :-80.4460  
##  1st Qu.: 14.00   1st Qu.: -50.16   1st Qu.: -0.3025  
##  Median : 30.00   Median :  17.66   Median : 16.5110  
##  Mean   : 54.84   Mean   :  13.29   Mean   : 16.4075  
##  3rd Qu.: 61.00   3rd Qu.:  50.01   3rd Qu.: 39.1067  
##  Max.   :155.00   Max.   : 179.22   Max.   : 78.8300

Collate countries and their names

Unfortunately, as soon as it comes to map the world, geopolitics come to the fore. This is reflected in the mismatches between country lists due to the different recognition statuses of these countries.

For one, there are distinctly more countries in the shp file:

length(world.spdf$NAME) 

## [1] 246

length(epi$country)

## [1] 179

These mismatches stem mainly from different naming schemes.

sort(setdiff(world.spdf$NAME, epi$country))

##  [1] "Åland Islands"                            
##  [2] "American Samoa"                           
##  [3] "Andorra"                                  
##  [4] "Anguilla"                                 
##  [5] "Antarctica"                               
##  [6] "Aruba"                                    
##  [7] "Bermuda"                                  
##  [8] "Bouvet Island"                            
##  [9] "British Indian Ocean Territory"           
## [10] "British Virgin Islands"                   
## [11] "Burma"                                    
## [12] "Cape Verde"                               
## [13] "Cayman Islands"                           
## [14] "Christmas Island"                         
## [15] "Cocos (Keeling) Islands"                  
## [16] "Congo"                                    
## [17] "Cook Islands"                             
## [18] "Democratic Republic of the Congo"         
## [19] "Falkland Islands (Malvinas)"              
## [20] "Faroe Islands"                            
## [21] "French Guiana"                            
## [22] "French Polynesia"                         
## [23] "French Southern and Antarctic Lands"      
## [24] "Gibraltar"                                
## [25] "Greenland"                                
## [26] "Guadeloupe"                               
## [27] "Guam"                                     
## [28] "Guernsey"                                 
## [29] "Heard Island and McDonald Islands"        
## [30] "Holy See (Vatican City)"                  
## [31] "Hong Kong"                                
## [32] "Iran (Islamic Republic of)"               
## [33] "Isle of Man"                              
## [34] "Jersey"                                   
## [35] "Korea, Democratic People's Republic of"   
## [36] "Korea, Republic of"                       
## [37] "Lao People's Democratic Republic"         
## [38] "Libyan Arab Jamahiriya"                   
## [39] "Liechtenstein"                            
## [40] "Macau"                                    
## [41] "Martinique"                               
## [42] "Mayotte"                                  
## [43] "Micronesia, Federated States of"          
## [44] "Monaco"                                   
## [45] "Montserrat"                               
## [46] "Nauru"                                    
## [47] "Netherlands Antilles"                     
## [48] "New Caledonia"                            
## [49] "Niue"                                     
## [50] "Norfolk Island"                           
## [51] "Northern Mariana Islands"                 
## [52] "Palau"                                    
## [53] "Palestine"                                
## [54] "Pitcairn Islands"                         
## [55] "Puerto Rico"                              
## [56] "Republic of Moldova"                      
## [57] "Reunion"                                  
## [58] "Saint Barthelemy"                         
## [59] "Saint Helena"                             
## [60] "Saint Kitts and Nevis"                    
## [61] "Saint Martin"                             
## [62] "Saint Pierre and Miquelon"                
## [63] "San Marino"                               
## [64] "Somalia"                                  
## [65] "South Georgia South Sandwich Islands"     
## [66] "Svalbard"                                 
## [67] "Swaziland"                                
## [68] "Syrian Arab Republic"                     
## [69] "Taiwan"                                   
## [70] "The former Yugoslav Republic of Macedonia"
## [71] "Tokelau"                                  
## [72] "Turks and Caicos Islands"                 
## [73] "Tuvalu"                                   
## [74] "United Republic of Tanzania"              
## [75] "United States"                            
## [76] "United States Minor Outlying Islands"     
## [77] "United States Virgin Islands"             
## [78] "Wallis and Futuna Islands"                
## [79] "Western Sahara"                           
## [80] "Yemen"

sort(setdiff(epi$country, world.spdf$NAME))

##  [1] "Cabo Verde"               "Dem. Rep. Congo"         
##  [3] "Eswatini"                 "Iran"                    
##  [5] "Laos"                     "Micronesia"              
##  [7] "Moldova"                  "Myanmar"                 
##  [9] "North Macedonia"          "Republic of Congo"       
## [11] "South Korea"              "Tanzania"                
## [13] "United States of America"

Now, collate the naming schemes between the epi data.frame and the world.spdf SpatialPolygonsDataFrame with the helper function collateCountryNames().

epi.world.spdf <- collateCountryNames(epi, world.spdf)
list2env(epi.world.spdf, .GlobalEnv)

## <environment: R_GlobalEnv>

rm(epi.world.spdf)

Lastly, control if all EPI countries are present in the world map…

sort(setdiff(epi$country, world.spdf$NAME))

## character(0)

… as the set difference is now equal to zero (i.e., all items in epi$country are present also in world.spdf$NAME).

Countries without an EPI are miniature states, islands or states currently in an armed conflict (some of the island belong to a nation; e.g., Greenland to Denkmark and Svalbard to Norway. For now, keep it as is and later, these islands could possibly be merged into their mainlands).

sort(setdiff(world.spdf$NAME, epi$country))

##  [1] "Åland Islands"                       
##  [2] "American Samoa"                      
##  [3] "Andorra"                             
##  [4] "Anguilla"                            
##  [5] "Antarctica"                          
##  [6] "Aruba"                               
##  [7] "Bermuda"                             
##  [8] "Bouvet Island"                       
##  [9] "British Indian Ocean Territory"      
## [10] "British Virgin Islands"              
## [11] "Cayman Islands"                      
## [12] "Christmas Island"                    
## [13] "Cocos (Keeling) Islands"             
## [14] "Cook Islands"                        
## [15] "Falkland Islands (Malvinas)"         
## [16] "Faroe Islands"                       
## [17] "French Guiana"                       
## [18] "French Polynesia"                    
## [19] "French Southern and Antarctic Lands" 
## [20] "Gibraltar"                           
## [21] "Greenland"                           
## [22] "Guadeloupe"                          
## [23] "Guam"                                
## [24] "Guernsey"                            
## [25] "Heard Island and McDonald Islands"   
## [26] "Holy See (Vatican City)"             
## [27] "Hong Kong"                           
## [28] "Isle of Man"                         
## [29] "Jersey"                              
## [30] "Libyan Arab Jamahiriya"              
## [31] "Liechtenstein"                       
## [32] "Macau"                               
## [33] "Martinique"                          
## [34] "Mayotte"                             
## [35] "Monaco"                              
## [36] "Montserrat"                          
## [37] "Nauru"                               
## [38] "Netherlands Antilles"                
## [39] "New Caledonia"                       
## [40] "Niue"                                
## [41] "Norfolk Island"                      
## [42] "North Korea"                         
## [43] "Northern Mariana Islands"            
## [44] "Palau"                               
## [45] "Palestine"                           
## [46] "Pitcairn Islands"                    
## [47] "Puerto Rico"                         
## [48] "Reunion"                             
## [49] "Saint Barthelemy"                    
## [50] "Saint Helena"                        
## [51] "Saint Kitts and Nevis"               
## [52] "Saint Martin"                        
## [53] "Saint Pierre and Miquelon"           
## [54] "San Marino"                          
## [55] "Somalia"                             
## [56] "South Georgia South Sandwich Islands"
## [57] "Svalbard"                            
## [58] "Syrian Arab Republic"                
## [59] "Taiwan"                              
## [60] "Tokelau"                             
## [61] "Turks and Caicos Islands"            
## [62] "Tuvalu"                              
## [63] "United States Minor Outlying Islands"
## [64] "United States Virgin Islands"        
## [65] "Wallis and Futuna Islands"           
## [66] "Western Sahara"                      
## [67] "Yemen"

Add EPI and regions to sp file

world.spdf <- merge(world.spdf,
                    epi[, (colnames(epi) %in% c("country", "EPI.new", "region",
                                                "X2019"))], 
                    by.x = "NAME", by.y = "country")
names(world.spdf)

##  [1] "NAME"      "FIPS"      "ISO2"      "ISO3"      "UN"        "AREA"     
##  [7] "POP2005"   "REGION"    "SUBREGION" "LON"       "LAT"       "EPI.new"  
## [13] "region"    "X2019"

Encode population numbers as millions

world.spdf@data$POP2005[which(world.spdf@data$POP2005 == 0)] <- NA
world.spdf@data$POP2005 <- as.integer(world.spdf@data$POP2005) / 1000000 %>%
  round(2)

Choropleth map with leaflet

Create a color palette for the map

The EPI color scheme:

epi.cols <- colorNumeric(palette = "viridis", domain = world.spdf@data$EPI.new,
                         na.color = "gray", reverse = TRUE)

The region color scheme:

region.cols <- colorFactor(palette = "plasma", domain = world.spdf@data$region,
                           na.color = "gray")

Prepare the text for tooltips

tooltips <- paste("Country:", world.spdf@data$NAME, "<br/>",
                  "Region:", world.spdf@data$region, "<br/r>",
                  # "Area:", world.spdf@data$AREA, "<br/>",
                  "Population:", round(world.spdf@data$POP2005, 2), "M",
                  "<br/>",
                  "EPI:", world.spdf@data$EPI.new) %>%
  lapply(htmltools::HTML)

Render map

cp.map <- leaflet(world.spdf) %>%
  addTiles() %>%
  setView(lat = 20, lng = 0, zoom = 2.5) %>%
  addPolygons(fillColor = ~ epi.cols(EPI.new), stroke = TRUE, group = "EPI",
              fillOpacity = 0.9, color = "white", weight = 0.5,
              label = tooltips,
              labelOptions = labelOptions(style = list("font-weight" = "normal",
                                                       padding = "3px 8px"),
                                          textsize = "13px",
                                          direction = "auto")) %>%
  addPolygons(fillColor = ~ region.cols(region), stroke = TRUE,
              group = "Region", fillOpacity = 0.5, color = "white",
              weight = 0.5, label = tooltips,
              labelOptions = labelOptions(style = list("font-weight" = "normal",
                                                       padding = "3px 8px"),
                                          textsize = "13px",
                                          direction = "auto")) %>%
  addLegend(pal = epi.cols, values = ~ EPI.new, opacity = 0.9,
            title = "EPI", position = "topleft", group = "EPI") %>%
  addLegend(pal = region.cols, values = ~ region, opacity = 0.5,
            title = "", position = "topleft", group = "Region") %>%
  addLayersControl(overlayGroups = c("EPI", "Region"),
                   options = layersControlOptions(collapsed = FALSE)) %>%
  hideGroup("Region")
cp.map

Export map to html

This takes some time, so only run it if needed.

# saveWidget(cp.map, file = "index.html")

Analysis

Color scale for plotting

cols.region <- c("#0D0887", "#5402A3", "#8B0AA5", "#B93289", "#DB5C68",
                 "#F48849", "#FEBC2A", "#F0F921")

National comparison

The 32 performance indicators are displayed by country in boxplots sorted by their median and colored by region.

ggplot(epi.long[epi.long$Type == "Indicator", ],
       aes(x = reorder(country, EPI.new.value, FUN = median, na.rm = TRUE),
           y = EPI.new.value, fill = region)) +
  geom_boxplot(outlier.shape = 21, outlier.size = 1, alpha = 0.5) +
  scale_fill_manual(values = cols.region) +
  theme(legend.position = "top", legend.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.text = element_text(size = 5)) +
  ylab("2020 EPIs") +
  xlab("")

## Warning: Removed 446 rows containing non-finite values (stat_boxplot).

# ggsave("EPI_nations.pdf", width = 20, height = 8.27, bg = "transparent")

Regional comparison

The ?? indicators are pooled by region and are visually displayed in boxplots sorted by their median and colored by region.

ggplot(epi.long[epi.long$Type == "EPI", ],
       aes(x = reorder(region, EPI.new.value, FUN = median, na.rm = TRUE),
           y = EPI.new.value, fill = region)) +
  geom_boxplot(outlier.shape = 21, outlier.size = 1, alpha = 0.5) +
  scale_fill_manual(values = cols.region) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  ylim(0, 100) +
  ylab("2020 EPI") +
  xlab("")

# ggsave("EPI_regions.pdf", width = 11.29, height = 6.27, bg = "transparent")

EPI vs GDP

ggplot(epi.long[epi.long$Type == "EPI", ],
       aes(x = log(X2019 / 1000000), y = EPI.new.value, fill = region)) +
  geom_smooth(aes(group = 1), color = "gray", lty = 2,
              method = "lm", formula = y ~ x, se = FALSE, show.legend = FALSE) +
  geom_point(alpha = 0.5, pch = 21) +
  scale_color_manual(values = cols.region) +
  theme(legend.position = "top", legend.title = element_blank()) +
  ylim(0, 100) +
  ylab("2020 EPI") +
  xlab(expression(paste(log[e], "(2019 GDP)", " [US$]")))

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

# ggsave("EPI_GDP.pdf", width = 11.29, height = 6.27, bg = "transparent")

Structural equation model

In a structural equation model …

Merge datasets

epi <- merge(epi,
             world.spdf[, (names(world.spdf) %in% c("ISO3", "AREA",
                                                    "POP2005", "LAT"))],
             by.x = "iso", by.y = "ISO3")
names(epi)

##  [1] "iso"     "country" "code"    "EPI.new" "HLT.new" "AIR.new" "PMD.new"
##  [8] "HAD.new" "OZD.new" "H2O.new" "USD.new" "UWD.new" "HMT.new" "PBD.new"
## [15] "WMG.new" "MSW.new" "ECO.new" "BDH.new" "TBN.new" "TBG.new" "MPA.new"
## [22] "PAR.new" "SHI.new" "SPI.new" "BHV.new" "ECS.new" "TCL.new" "GRL.new"
## [29] "WTL.new" "FSH.new" "FSS.new" "RMS.new" "FGT.new" "CCH.new" "CDA.new"
## [36] "CHA.new" "FGA.new" "NDA.new" "BCA.new" "LCB.new" "GIB.new" "GHP.new"
## [43] "APE.new" "SDA.new" "NXA.new" "AGR.new" "SNM.new" "WRS.new" "WWT.new"
## [50] "region"  "X2019"   "AREA"    "POP2005" "LAT"

str(epi)

## 'data.frame':    179 obs. of  54 variables:
##  $ iso    : chr  "AFG" "AGO" "ALB" "ARE" ...
##  $ country: chr  "Afghanistan" "Angola" "Albania" "United Arab Emirates" ...
##  $ code   : int  4 24 8 784 32 51 28 36 40 31 ...
##  $ EPI.new: num  25.5 29.7 49 55.6 52.2 52.3 48.5 74.9 79.6 46.5 ...
##  $ HLT.new: num  20 20.4 44.5 55.2 60.2 43.5 55.5 91.6 88 32.7 ...
##  $ AIR.new: num  17.7 26.8 41.2 48.6 56.9 36.3 57.5 98.2 81.3 24.9 ...
##  $ PMD.new: num  25.4 32.6 44.2 13.9 54.7 22.1 44.6 100 73.9 12.7 ...
##  $ HAD.new: num  7.3 17.8 34.6 100 60.9 57 69.8 98.8 96.4 40.7 ...
##  $ OZD.new: num  16.9 34.4 60.1 18.8 48.7 26.1 100 73.8 42.2 33.1 ...
##  $ H2O.new: num  28 12.8 54 67.2 64.7 57.2 50 87 94.7 45.5 ...
##  $ USD.new: num  28.5 14.1 59.7 90.5 72.7 50.7 57.5 96 86.7 44.3 ...
##  $ UWD.new: num  27.7 11.9 50.2 51.6 59.3 61.5 45 81 100 46.4 ...
##  $ HMT.new: num  0 37.3 46.1 54.3 73.1 50.8 60.5 77.2 91.7 41 ...
##  $ PBD.new: num  0 37.3 46.1 54.3 73.1 50.8 60.5 77.2 91.7 41 ...
##  $ WMG.new: num  0 0 0 26.8 44.6 0 75.1 77.3 97.2 0 ...
##  $ MSW.new: num  0 0 0 26.8 44.6 0 75.1 77.3 97.2 0 ...
##  $ ECO.new: num  29.2 35.9 52 55.9 46.8 58.1 43.8 63.8 74 55.7 ...
##  $ BDH.new: num  21.9 39.3 68.2 80.9 49.1 79.2 55.7 83.7 85.5 56.9 ...
##  $ TBN.new: num  0.6 34.1 100 90.4 41.7 100 90.3 89 100 50.7 ...
##  $ TBG.new: num  1.3 36.6 100 90.4 56.7 100 87.3 82.5 100 62.2 ...
##  $ MPA.new: num  NA 0 14.7 100 43 NA 1.2 100 NA NA ...
##  $ PAR.new: num  5 29.8 38.6 8.9 16.2 19 53.8 40.4 63.7 11.7 ...
##  $ SHI.new: num  90.7 88.3 82.1 81.3 69.4 92.4 NA 87.3 84.7 96 ...
##  $ SPI.new: num  18 66.6 92.2 80.3 64.3 74.5 NA 92.4 81.2 77.3 ...
##  $ BHV.new: num  57.7 67.2 39.6 76.8 58.9 48 34 73.4 54.7 44.4 ...
##  $ ECS.new: num  93.6 33.7 43.2 100 29.4 81.5 38.8 27.9 35.6 80.5 ...
##  $ TCL.new: num  95 31.2 41.1 NA 27.6 82.4 35.4 24.5 32.5 80.2 ...
##  $ GRL.new: num  60.9 49.3 70.2 100 47.3 45.8 39.8 44 69.9 65.2 ...
##  $ WTL.new: num  100 63.5 54.4 100 44.5 100 100 72.9 57.3 100 ...
##  $ FSH.new: num  NA 14.7 8.9 14.5 3.2 NA 5 3.4 NA NA ...
##  $ FSS.new: num  NA 16.9 NA 13.5 0.7 NA 1.4 0.7 NA NA ...
##  $ RMS.new: num  NA 21.3 16.2 15.6 6.2 NA 8.6 NA NA NA ...
##  $ FGT.new: num  NA 4.4 0.4 NA 2.6 NA NA 6.5 NA NA ...
##  $ CCH.new: num  22.2 49 56.8 38.9 60.2 46.7 58.5 70.4 71.3 48.6 ...
##  $ CDA.new: num  0 30.3 41.2 26.5 48.2 39.3 48.4 66.7 62.8 42.4 ...
##  $ CHA.new: num  66.1 100 79.2 52.7 88.1 54.7 100 85.8 100 43.8 ...
##  $ FGA.new: num  NA 89.3 89.4 92.6 90.7 89.4 89.8 90.9 87.8 88.9 ...
##  $ NDA.new: num  34.1 63.2 61.9 41.7 58 0 77.8 70.1 66.9 60.9 ...
##  $ BCA.new: num  15.7 7.1 100 74.4 90.3 57.6 0 56.9 100 37.4 ...
##  $ LCB.new: num  100 27.2 78.4 15 57.8 97.3 79.7 89.3 76.7 100 ...
##  $ GIB.new: num  50.6 53.9 36.4 19.8 36.2 23 22 62.5 42.7 27.7 ...
##  $ GHP.new: num  100 60.3 58.3 0 29.8 60.3 35.2 0 23 40.3 ...
##  $ APE.new: num  0 2.6 100 83.4 59.3 46.8 36.6 90.3 100 100 ...
##  $ SDA.new: num  0 5.2 100 100 67.3 50.6 28.2 92.4 100 100 ...
##  $ NXA.new: num  0 0 100 66.7 51.3 43 45 88.2 100 100 ...
##  $ AGR.new: num  51 29.3 37.6 13.7 78.4 57.6 5.1 49.2 68 64.3 ...
##  $ SNM.new: num  51 29.3 37.6 13.7 78.4 57.6 5.1 49.2 68 64.3 ...
##  $ WRS.new: num  0 0 2.7 76.8 5.9 8.8 1.3 92.7 94 3.8 ...
##  $ WWT.new: num  0 0 2.7 76.8 5.9 8.8 1.3 92.7 94 3.8 ...
##  $ region : chr  "Southern Asia" "Sub-Saharan Africa" "Eastern Europe" "Greater Middle East" ...
##  $ X2019  : num  1.93e+10 8.88e+10 1.53e+10 4.21e+11 4.45e+11 ...
##  $ AREA   : int  65209 124670 2740 8360 273669 2820 44 768230 8245 8260 ...
##  $ POP2005: num  25.07 16.1 3.15 4.1 38.75 ...
##  $ LAT    : num  33.7 -12.3 41.1 23.5 -35.4 ...

Scale the variables

epi[, c(4, 51:54)] <- apply(epi[, c(4, 51:54)], 2, scale)

Specify the model

epi.gdp <- 'EPI.new ~ X2019 + POP2005 + LAT
  X2019 ~ LAT + POP2005
  POP2005 ~ AREA'

Run the model and inspect the results

fit.epi.gdp <- sem(epi.gdp, data = epi)
summary(fit.epi.gdp, rsq = TRUE)

## lavaan 0.6-7 ended normally after 13 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of free parameters                          9
##                                                       
##                                                   Used       Total
##   Number of observations                           172         179
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                26.066
##   Degrees of freedom                                 3
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   EPI.new ~                                           
##     X2019             0.292    0.078    3.733    0.000
##     POP2005          -0.271    0.076   -3.560    0.000
##     LAT               0.509    0.062    8.245    0.000
##   X2019 ~                                             
##     LAT               0.120    0.059    2.017    0.044
##     POP2005           0.582    0.059    9.830    0.000
##   POP2005 ~                                           
##     AREA              0.459    0.068    6.775    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .EPI.new           0.657    0.071    9.274    0.000
##    .X2019             0.622    0.067    9.274    0.000
##    .POP2005           0.815    0.088    9.274    0.000
## 
## R-Square:
##                    Estimate
##     EPI.new           0.357
##     X2019             0.372
##     POP2005           0.211

Inspect modification indices

The modification indices suggest–among others–that there is an unspecified cause that drives GDP and population numbers. In SEM, this unspecified result is modeled by a correlation and encoded with this double tilde ~~.

modificationindices(fit.epi.gdp, minimum.value = 3)

##        lhs op     rhs     mi    epc sepc.lv sepc.all sepc.nox
## 15   X2019 ~~ POP2005 24.039 -0.581  -0.581   -0.816   -0.816
## 18   X2019  ~    AREA 24.039  0.327   0.327    0.334    0.329
## 20 POP2005  ~   X2019 20.315 -0.814  -0.814   -0.797   -0.797
## 27    AREA  ~   X2019 21.889  0.494   0.494    0.484    0.484

fit.epi.gdp.up <- update(fit.epi.gdp, add = "X2019 ~~ POP2005")
summary(fit.epi.gdp.up, rsq = TRUE)

## lavaan 0.6-7 ended normally after 23 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of free parameters                         10
##                                                       
##                                                   Used       Total
##   Number of observations                           172         179
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                 0.211
##   Degrees of freedom                                 2
##   P-value (Chi-square)                           0.900
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   EPI.new ~                                           
##     X2019             0.292    0.078    3.732    0.000
##     POP2005          -0.271    0.076   -3.560    0.000
##     LAT               0.509    0.062    8.238    0.000
##   X2019 ~                                             
##     LAT               0.109    0.055    1.981    0.048
##     POP2005           1.144    0.159    7.181    0.000
##   POP2005 ~                                           
##     AREA              0.459    0.068    6.775    0.000
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##  .X2019 ~~                                            
##    .POP2005          -0.580    0.153   -3.799    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .EPI.new           0.657    0.071    9.274    0.000
##    .X2019             0.948    0.211    4.488    0.000
##    .POP2005           0.815    0.088    9.274    0.000
## 
## R-Square:
##                    Estimate
##     EPI.new           0.358
##     X2019             0.044
##     POP2005           0.211

Visualize the model

lavaanPlot(model = fit.epi.gdp.up,
           node_options = list(shape = "box", color = "gray",
                               fontname = "Helvetica"),
           edge_options = list(color = "black"),
           coefs = TRUE, covs = TRUE, stars = c("covs", "regress"),
           labels = list(AREA = "Area", POP2005 = "Population", X2019 = "GDP",
                         LAT = "Latitude", EPInew = "EPI"))